{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "power_method (generic function with 1 method)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function power_method(A,d;T=5000, tol=big(1e-100))\n",
    "    (n,n) = size(A)\n",
    "    x_all = big.(qr(randn(n,d))[1])\n",
    "    _eig  = zeros(BigFloat, d)\n",
    "    for j=1:d\n",
    "        x = x_all[:,j]\n",
    "        x /= norm(x)\n",
    "        for t=1:T\n",
    "            x = A*x\n",
    "            if j > 1\n",
    "                x -= sum(x_all[:,1:(j-1)]*diagm(vec(x'x_all[:,1:(j-1)])),2)\n",
    "            end\n",
    "            nx = norm(x)\n",
    "            x /= nx\n",
    "            if abs(nx - _eig[j]) < tol\n",
    "                println(\"\\t Done with eigenvalue $(j) at iteration $(t) at $(Float64(abs(nx - _eig[j]))) \")\n",
    "                break\n",
    "            end\n",
    "            _eig[j]    = nx\n",
    "        end\n",
    "        x_all[:,j] = x \n",
    "    end\n",
    "    return (_eig, x_all)\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 200.773\n",
       " 180.748\n",
       " 160.172"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A0 = randn(50,50)\n",
    "A  = A0'A0\n",
    "sort(eig(A)[1], rev=true)[1:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\t Done with eigenvalue 1 at iteration 850 at 0.0 \n",
      "\t Done with eigenvalue 2 at iteration 717 at 0.0 \n",
      "\t Done with eigenvalue 3 at iteration 953 at 0.0 \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " 200.773\n",
       " 180.748\n",
       " 160.172"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(eig2, v) = power_method(A,3)\n",
    "Float64.(eig2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([0.0016145, 0.0337414, 0.281478, 0.481398, 0.868266, 1.33826, 2.30098, 2.97151, 3.55055, 4.85673  …  95.3252, 101.001, 117.442, 119.162, 122.824, 129.312, 136.324, 145.041, 167.855, 191.365], [0.142506 0.0446259 … 0.0574164 -0.0944285; 0.122999 0.00241461 … -0.0799893 0.260021; … ; -0.137663 -0.0717979 … -0.0859499 -0.0136689; -0.0530304 0.187545 … 0.129448 -0.182961])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eig(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
